% version 1.71
% date 2022/1/4 3:16 pm
% auther gxliu

function FplotSeiq(yT)

	% number of days
	rb = 30;

	% I_warn: if I >= Iwarn ,then goverment will take action.
	I_warn = 100;

	% j = a2 a3 b1 b2, public
	j = [0.8756, 0.8315, 0.1000, 0.1054];

	%% First


	% first stage domin
	tspan1 = [0 rb];

	% r a1
	ra1 = [2, 1/1.51];
	% y0 = S, E, I, Q
	y01 = [10000000, 1, 0, 0];

	% first stage
	function dydt1 = vdp46(t, y, j)

		r = 20;
		a1 = 1/5.1;

		% a2 = 0;
		% a3 = 0;
		b1 = j(3);
		% b2 = j(4);

		S = y(1);
		E = y(2);
		I = y(3);
		Q = y(4);

		dydt1 = [ -r * b1 * E;
			   r * b1 * E - E * a1;
			   a1 * E;
				0];
	end

	[t1, y1] = ode45(@(t,y) vdp46(t,y,j),tspan1, y01);
	% show t1
	% t1 

	sy1 = size(y1);

	% anst, between f-stage to s-stage
	anst = -1;
	for ntimes = 1:length(t1)
		if (y1(ntimes,3) >= I_warn && anst == -1)
			anst = ntimes;
		end 
	end 

	% show anst and Break_time
	anst 
	% Break_time is a turning point when the government finds Iwarn cases(maybe one or more) and starts taking action
	Break_time = t1(anst,1)

	plot(t1(1:anst,:), y1(1:anst,2), t1(1:anst,:), y1(1:anst,3), t1(1:anst,:), y1(1:anst, 4));
	legend("E", "I", "Q");

	% second stage

	tspan2 = [Break_time rb];

	y02 = [y1(anst, 1), y1(anst, 2), y1(anst, 3), y1(anst, 4)];

	% show y02
	% y02(1) % 1e7
	% y02(2) % 11
	% y02(3) % 1.0992
	% y02(4) % 0

	function dydt2 = vdp47(t, y, j, ra1)

		r = ra1(1);
		a1 = ra1(2);

		a2 = j(1);
		a3 = j(2);
		b1 = j(3);
		b2 = j(4);

		S = y(1);
		E = y(2);
		I = y(3);
		Q = y(4);

		dydt2 = [-r*( b1*E + b2*I );
			  r*( b1*E + b2*I ) - E*( a1+a2 );
			  a1*E - a2*I;
			  a3*E + a2*I];
	end 

	[t2, y2] = ode45(@(t,y) vdp47(t,y,j,ra1),tspan2,y02);
	% show t2
	% t2
	
	t3 = cat(1, t1(1:anst,1), t2);
	st3 = size(t3);
	t3;

	y3 = cat(1, y1(1:anst,:), y2);
	sy3 = size(y3);

	syT = size(yT);
	tyT = [1:syT(1)];

	% plot(t3, y3(:,2), t3,y3(:,3), t3,y3(:,4));
	% legend("E","I","Q");
	plot(t3, y3(:,2), t3,y3(:,3), t3,y3(:,4), tyT,yT,'o');
	legend("E","I","Q","yT");

	% plot(t3, y3(:,4));
	% plot(t3, y3(:,4));

end
